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ABSTRACT. This paper describes a number of ways to improve on the standard 
method for measuring the two-point correlation function of large scale structure in the 
Universe. Issues addressed are: (1) The problem of the mean density, and how to 
solve it; (2) How to estimate the uncertainty in a measured correlation function; (3) 
Minimum variance pair weighting; (4) Unbiased estimation of the selection function 
when magnitudes are discrete; (5) Analytic computation of angular integrals in 
background pair counts. 


1. The Mean Density Problem 

It is widely thought that the accuracy of the correlation function § is fundamentally 
limited by uncertainty in the mean density. Actual, this notion is false (§1.2), although 
it is true for the commonly used estimator of t, (§1.1). 

l.l The Problem 


The statistic commonly used to estimate % from a catalog of galaxies is (e.g. Davis and 
Peebles 1983) 


= <MVxW> 
^ est <NW><N> 


( 1 ) 


where N represents real galaxies, W is the catalog window, and <> denotes averaging 
over all points in the catalog; for <NN> and <NW> the averaging is over all pairs 
lying in an interval of separations r . The observed galaxy density N is the true galaxy 
density n times the catalog window W (strictly, observed discrete galaxies are taken to 
be a Poisson process superimposed on this). In terms of the true galaxy overdensity 
8 = (n -n)/n, where n is the true mean galaxy density of the Universe, the observed 
galaxy density N is 

N=nW(l+8) . (2) 


To see how good the estimate (1) of_£ is in a realistically unfair sample, 
introduce the following notations (3)-(5). Let 8 be the mean overdensity in the catalog 


<1Y<5> 
o = 

<w> 


( 3 ) 
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( 4 ) 


and let y denote the galaxy-catalog correlation function 

_ <W j5j W 2 > 

¥n = <w x w 2 > 


In a fair sample, the mean overdensity 5 and the galaxy-catalog correlation function y 
would be identically zero, but they are not necessarily zero in reality. Let | denote the 
windowed galaxy-galaxy correlation function 


<W 1 8 1 W 2 5 2> 
<W 2 > 


(5) 


While | is not necessarily equal to the true correlation function £, of the Universe, it is 
at least the ‘true’ correlation function of the sample, which is presumably the next best 
thing. 


In terms of 8, y , and |, the standard estimate (1) of the correlation function £, is 


SestO") = 


_ $(r) + y(r) - 8 - y(r)8 


(1 + y(r ))(1 + 5) 


( 6 ) 


The problem with equation (6) is that it contains not-necessarily-vanishing terms 
y(r)-8 which are of first prder in overdensity 5, whereas the thing you want, the 
sample correlation function is of second order in 8. This is a severe drawback of 
the estimator (1) for l; in the linear regime of small 5. 


1.2 A Solution, part l 


A better estimate of £ is 

„ . <NNxWW> 

Sest “ 2 

<NW> 1 


(7) 


in which the brackets <> in both numerator and denominator denote averaging over 
pairs in an interval of separations r. In terms of the galaxy-catalog correlation 
function y and the sample galaxy-galaxy correlation function | defined by equations 
(4) & (5), the estimate (7) is 


Z M ir)= 

(1 + y(r)) 2 


( 8 ) 


which differs from the sample correlation function £ only by terms which are of 
second order in overdensity 8. 

The advantages of the estimator (7) over the standard estimator (1) for £ are: 


(a) Accuracy, especially in the large scale, linear regime; 

(b) Reliability in the presence of unfairness, especially with not-unbiased (e.g. 
minimum variance, §3) pair weightings; 

(c) Peace of mind: there is no need to measure the mean density <N>/<W> as a 
separate operation; equation (7) specifies that the ‘correct’ mean density to use in 
place of the <N>/<W> in equation (1) is <NW> 2 /<WW>, a quantity which it is 
to be noted varies with separation r . 


16 


1.3 SOLUTION, PART 2 


The y/ 2 term in equation (8) represents large scale variance which is inevitably missing 
in a finite catalog; its presence is symptomatic of the familiar problem that using the 
sample mean leads to an underestimate of the sample variance. Although the galaxy- 
catalog correlation y/ should be zero in the mean over many samples, its variance 
<yrh> should be positive in the mean. Physically, <y/ 2 > represents the mean fractional 
excess of galaxies clustered around a galaxy on the scale of the catalog. If one 
imagines evaluating the correlation function by sitting on galaxies and counting 
neighbors, the mean density at infinity should be determined not from the total number 
of galaxies in the catalog volume, but rather from the number of galaxies less the 
mean excess of galaxies clustered around a galaxy. 

One way to correct for the missing variance is to use an alternate estimator 


£est 


<NNxWW> 


- 1 


(9) 


<AW> 2 - <A(AW) 2 > 
which in terms of f and yr and its variance <yr 2 > is 

U(D = ■ (10) 

(1 + y/(r)) 2 - <yr(r) 2 > 

The variance <A(AW) 2 > in equation (9) may be computed from the fluctuations in 
<AW>, using methods similar to those described in §2. Another way to correct for 
missing large scale variance is given by Hamilton (1993). 


2. Estimating Errors in the Correlation Function 


2.1 mathematics 


The window W can be imagined as a set of weights W t attached to every tiny volume 
element of the Universe. For an observed subsample, the weights W i are nonzero only 
over the observed region. For the entire population, the Universe, the weights VF ( - ipop 
are finite everywhere, but infinitesimal compared to the sample weights W { . 

The correlation function £(VF,) measured in a sample differs from the true 
correlation function ^(V^ J ip0 p) by an error 


A£ = W)- . 

Expanding the error as a Taylor series to second order in the weights gives 


( 11 ) 


AS=2(^-^,pop) 






a 2 s 


pop 


+ yZW -"Wa w#Wj 


.(12) 


Ipop 


Using the facts that (a) £(VF,) is unchanged by rescaling (b) ^ is a 

quadratic function of the weights W i (at least for the estimator [7]), and (c) VP t - pop is 
infinitesimal compared to W it eliminates most of the terms in equation .(12), reducing 
it to 


Af = — T WW - — — — 
? 2 1 J dW s dWj 


pop 


which then reduces further to 


(13) 


17 


A£= X 

distinct ij 


w, 




" Wij 


[pop 


(14) 


where = WiWj. Expression (14) makes clear the fact that the error is truly a 
derivative with respect to pairs. Approximating the population derivative of £ in (14) 
by the sample derivative, and again using the fact that <^(W ( ) is quadratic in W { , 
permits equation (14) to be rearranged as a sum over volume elements i rather than 
pairs ij : 


= I AS,- 


with 


AS; = — Wi 

2 1 dWi 


(15) 


Note the important factor of V 2 in equation (15), which in effect causes pairs to be 
counted once, not twice. The variance of S is then 

<Af> = £ A?, -A £,j . (16) 

ij 

Characteristically, the variance (16) increases as pairs ij of greater and greater 
separation are included, reaches a maximum, then declines to exactly zero when all 
pairs are included. The declining to zero is a consequence of approximating the 
population derivatives of S with the sample derivatives of S- To solve the problem, 
only pairs ij separated by some finite distance should be included. 

For the estimator (7), the contribution AS,- to the error in S from volume element 
i is 


A = (1 + §„.) 


<N { N> 

<NN> 


<N i W> <AW; > <W t W> 
<NW> ~ <NW> + <WW> 


(17) 


2.2 Suggested step-by-Step Error analysis 

(a) Divide the catalog into many subregions / . 

(b) Estimate £, from equation (7), and compute A for each subregion from equation 
(17). 

(c) Compute the variance <A£ 2 > from equation (16), including pairs ij of subregions 
of greater and greater separation, until the variance reaches a maximum. 

(d) The 1 -sigma error in <* is the square root of this variance. 


3. Minimum Variance Pair Weighting 
3.1 Mathematics 


An unbiased estimate of the correlation function % is gotten in principle by weighting 
each point inversely with the selection function O at the point, so that all volume 
elements count equally. Unfortunately this leads to a noisy estimate of % from regions 
where the selection function O is small. To the extent that the selection function is 
uncorrelated with the true galaxy density, the most accurate estimate of t, is obtained 
by reweighting the unbiased weighting of pairs inversely with the variance <A^ 2 > of £, 
so that a (real or background) pair 12 is weighted 


w 12 = 


1 

0j0 2 <A^ 2 > 


(18) 
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If the only source of uncertainty in % comes from the fact that the sample is a finite 
subset of the Universe, then the expected covariance between <Jj’s at separations r 12 
and r 34 is 


<A«^j2 A^34> — 1 <5 2 5354> <5i<52><53(54> 

^ J (^13^24 + ^14^23 + H 1234) ^3 (19) 


the integral being carried over all possible separations of point 3 from point 1. An 
important point to notice is that the £’s in the integrand of (19) have delta functions at 
zero separation, because of the discreteness of galaxies. These delta functions cause 
equation (19) to take the general form 

<A£ 2 > ~ O -2 + 20 _1 7 + K (20) 

in a region where the selection function is O. The O -2 term in equation (20) comes 
the case where pairs 12 and 34 coincide, the O -1 term from cases where 12 and 34 
share a point in common, and the constant term from cases where 12 and 34 are 
disjoint pairs. Equation (20) yields the pair weighting 


w ]2 = 


1 

1 + 207 + /CO 2 


( 21 ) 


Generally one is interested not in % at some precise separation, but rather 
averaged over some range of separations; or one might be interested in the power 
spectrum, or the harmonics of £, or such like. In that case equation (19) should be 
integrated with the desired kernel functions over the desired ranges of separations. 
The result is again equations of the form (20) and (21), but with different values of the 
coefficients 7 and K . The bad news is that a calculation of 7 and K from equation 
(19) is generally tricky and uncertain, and in any case suspect because an 
observationally measured £, may be subject to other sources of uncertainty ignored in 
equation (19). The good news is, the calculation is unnecessary. A practical solution 
is to proceed empirically, using the form (21) as a template for an approximate pair 
weighting, the free parameters of which would be determined empirically by 
minimizing the observed variance <A£ 2 > computed for example using the method of 
§2. This is the approach suggested in §3.2 below. 


3.2 Suggested Near minimum Variance Weighting 


A simple approximation to the minimum variance pair weighting which should not be 
too bad in practice would be to weight every (real and background) pair 12 by 
W 12 = H T W 2> the weight w,- at each point i where the selection function is O, being 


Wi = 


1 

1 + 0,7 


( 22 ) 


The quantity 7 in equation (22) is likely to be different for different pair separations r . 
Consideration of the behavior of the integral (19) suggests that a reasonable guess 
would be to take 


7 = Cr c (23) 

with c = 3-y if ^ x r~ r . The free parameters C and c in (23) would be determined 
empirically by varying them until the computed variance <A£ 2 > of or of whatever 
integral over £, is the quantity of interest, is minimized. Clearly the approximation 
(23), and perhaps also (22), could be refined if deemed necessary. 


19 


4. The Selection Function and the Discreteness of Magnitudes 

4.1 PROBLEM 

Turner’s (1979) classic inhomogeneity-insensitive method of measuring the selection 
function is found empirically to be sensitive to bin size. 

4.2 DIAGNOSIS 

Strauss, Yahil & Davis (1991) correctly attribute the sensitivity of Turner’s method to 
bin size to the fact that magnitudes given in catalogs are discretely, not continuously, 
distributed. 

The periodicity of listed magnitudes (e.g. Zwicky gives magnitudes mostly to 0.1) 
translates into a periodic ripple in the derived selection function (with period 0.1 
magnitudes in Zwicky’s case). 

4.3 SOLUTION 

The fix is simple: just choose a bin size equal to the period (0.1 magnitudes in 
Zwicky’s case), so the periodic ripple has no effect. Note that the worst possible 
binning is exactly half a period, since successive samplings of the selection function 
are then maximally out of phase. 


5. Better Backgrounds 

Oftentimes a catalog window is the product of a radial selection function and an 
angular window which is one inside, zero outside, a boundary composed of a set of arc 
segments. Angular integrals in such a case can be done analytically, and backgrounds 
can then be integrated quickly and accurately. 

[At the conference, a binder was exhibited containing listings of Fortran code to 
compute angular integrals analytically. The relevant mathematics is given by Hamilton 
(1993, Appendix).] 
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